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ABSTRACT 

Here we introduce CRASH3, the latest release of the 3D radiative transfer code CRASH. 
In its current implementation CRASH3 integrates into the reference algorithm the code 
Cloudy to evaluate the ionisation states of metals, self-consistently with the radiative 
transfer through H and He. The feedback of the heavy elements on the calculation of 
the gas temperature is also taken into account, making of CRASH3 the first 3D code for 
cosmological applications which treats self-consistently the radiative transfer through 
an inhomogeneous distribution of metal enriched gas with an arbitrary number of 
point sources and/or a background radiation. The code has been tested in idealized 
configurations, as well as in a more realistic case of multiple sources embedded in a 
polluted cosmic web. Through these validation tests the new method has been proven 
to be numerically stable and convergent. We have studied the dependence of the results 
on a number of physical quantities such as the source characteristics (spectral range 
and shape, intensity), the metal composition, the gas number density and metallicity. 
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1 INTRODUCTION 

During the last decade, observational and theoretical stud- 
ies constraining the nature of the intergalactic medium 
(IGM) have shown that metals are a pervasive component 
of the baryonic budget of our universe and that they are 
associated with a wide range of hydrogen column density 
(Nm) systems typically identified in quasars (QSO) absorp- 
tion s pectra at different redshift dMever and Yorklll987l; iLul 



n sp 

199l l; ISongaila and Cowid 119961; Cowie and Songaila 



Ellison et al 



Aracil et aL, 
Meiksinl 



2003 



200d: ISongailall200ll; jschave et alj|200d. 
2004 |Pieri et al.l | 2004 iBecker et al.ll2011l ; <ils<> 
2009 for a review) 



1998 



High Hi column density absorbers, with Nm ^ 2 ■ 
10 20 cm~ 2 , are easily identified in these spectra for the 
presence of strong damping wings, and are classified as 
Damped Lya systems (DLAs). Column densities in the 
range 1.6 • 10 17 cm" 2 Nm ^ 10 20 cm -2 are instead as- 
sociated with Lyman Limit Systems (LLSs) and spectro- 
scopically identified by their strongly saturated Lya lines. 
Both systems show C iv lines as well as many ions in lower 
ionisation states (Mnn, Sin, Fell), and are typicall y asso- 



The presence of metals in LLSs and DLAs can be in- 
terpreted as a natural product of the stellar nucleosynthesis 
acting therein. LLSs are in fact identified as clouds in the 
galactic halos, while high redshift (z ~ 3) DLAs are believed 
to be the progenitors of the present-day galaxies. 

Advances in high resolution spectroscopy re- 
vealed weak metal absorption lines in the Lya forest 
(Nm ~ 10" 17 cm" 2 ), with Civ detected in most of 
the systems with Nm J? 10 15 cm~ 2 and in more than 
half of the systems with Nm > 10 14 cm" 2 (|Tvtler et all 
1 19951 ; ISongaila and Cowid Il996h . The typical metal- 
licities of the Lya forest are estimated in the range 
10" 4 Z© < Z < 10" 2 Z© (|Simcoe et alj |2004 ). The sub- 
sequent discovery of a metal li c component in less dense 



regions (ICowie and Songailal 19981; Ellison et al. l200d; 



dated with metallicities of Z ~ 10 2 Zq (LLS s; e.g. Steide] 



1990h and 10~ 2 Zp) < Z < 0.3Z W (DLAs; e.g. ISteide 



Hellsten et al.lll997l; ITiauch et al.1ll997l ; iPettini et aL 
Kulkarni et alj|2005l h 



1990 



1999 
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I Schave et all l200d . [20031 ; lAracil et all 12004 iPieri et all 
2004 ) has been interpreted as the evidence of efficient 
spreading mechanisms which are able to transport the 
metals far away from their production sites and pollute low 
density regions. 

The redshift evolution of the gas metallicity has been 
extensively investigated. In the redshift range 1.5 < z < 4, 
C iv and Si iv doublets are the main tracers of the IGM 
metallicity, because their rest frame wavelength is larger 
than the Lya and the lines cannot be confused with those 
of the forest. In this range the colu mn density dis tribu- 
tion of Civ seems to remain constant (|Songailal l2001'). Al- 
though a consensus has not been reached yet, a decline in 
the abundance of Civ above z ~ 4.5 is reported by different 
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groups dRvan- Weber et~al1 120091 : iBecker et ail 120091 . l201ll : 
ISimcoe et all l201lh . TheSiiv doublet AA1393. 76, 1402.77 
A has been investigated in QSO li ne su rveys in the red- 
shift range 1.5 < z < 5.5. ISongailal (|200lh measured fisilV) 
the Si iv mass density relative to the critical density, in the 
redshift range 2 < z < 5.5, rinding a roughly constant 
value in the range 2 < z < 4.5 for absorbers with col- 
umn densities 10 13 cm~ 2 ^ iVsnv ^ 10 15 cm -2 , while fisirv 
may have increased by an order of magnitude in 4.5 < 
z < 5.5. These abu ndance trends have been confirmed by 
subsequent studies ( Boksenberg et al.l 120031 ; ISongailal 120051 : 
IScannapieco et alj|2006h . 

At higher redshifts, a critical diagnostic role could 
be played instead by ions with lower ionisat ion states 
like C 11, Si 11 and O n. probing colder gas (|Ohl |2002| ; 
iFurlanetto fc Loeblliool ; IBecker et al.ll201lft . 

Below z~la larger sample of data is available also for 
the oxygen component which is the major tracer of the IGM 
metallicity. The intere sted reader can find m ore information 
in the recent study of ICooksev et all (|201ll ). and references 
therein. 

Beyond this general consensus, many observational data 
remain controversial, as well as their interpretation which is 
of primary importance, e.g. in constraining different enrich- 
ment scenarios. The metal abundance, the number of ioni- 
sation states and the distributi on in space and time are still 
subjects of intense debate (see |Petitiean|[200ll . 120051 ). 

The theoretical interpretation of the data is of par- 
ticular relevance because observations of metals in dense 
regions, where stellar nucleosynthesis is active, provide a 
record of the star formation history, while observations in 
the IGM can be a good indicator of the galactic winds 
efficiency, the velocity structure of t he IGM and more in 
general of the enrichment mechanism { Gncdin and O strikei 
1997l:lFerrara et al.ll2000l;ICen and Bryan 2001; Madau et al 



20011 ). In fact, different scenarios have been considered in 
the literature, as an early enr i chmen t by the first gen- 
eration of stars (Mad a~u"et al.l l200lft . a continuous en- 



richment dGnedin and Ostriker! Il997i : berrara et"afl l200d : 
IScannapieco et al.l |2002| ) . or a late enric hment coinciding 
with the star formation peak at z ~ 2 — 4 (|Adelberger et ail 
2002). In addition, the determination of metal abundances 
could also constrain the effic i ency of the gas cooling function 
^Sutherland fc Dopital Il993[ : iMaio et alj 120071 ; ISmith et al.l 
120081 : Wiersma et alj|2009aft and the formation of massive 
galaxies IjThacker et al.l l2002). 

For these reasons, several theoretical schemes for 
metal production and spreading have been developed, 
based on bo th semi-analytic models (e.g. IScannapieco et aF 
20021 . 120061) and numerical simulations llMosconi et al 



200 



Scannapieco ct al. 2005; Oppenheimer & Dave 2006; 
Dubois fc Tevssierll2007l: IWiersma et al.l l2009bl: IMaio et al 
2010| : ISchave et alj |20ld: IShen et al.ll20ld; iTornatore et ~ 



201C : IWiersma et aij|20ld ; IMaio et al.ll201ll ; I Wiersma et a! 



20111 ). In the most advanced approaches the ionisation state 



of the metals is calculated assuming the presence of a uni- 
form UV background, which is then used as energetic in- 
put for photo-ionisation codes computing the metal ionisa- 
tion states at the e quilib rium (I Oppenheimer fc Pavel 120061 ; 
lOppenheimer et al. 2009. 2012ft . The results of this approach 
are affected by the uncertainties associated to the assump- 
tions on the shape and intensity of the radiation field, which 



is not calculated self-consistently from the radiative trans- 
fer across the inhomogeneous gas distribution. Many studies 
suggest in fact that shadowing, filtering and self-shielding 
induce deviations in the shape and intensity of the back- 
ground with respect to models in which the effects of the 
radiative transfer are neglected (|Maselli fc Ferrarall2005l and 
references therein). 

Fluctuations in the photo-ionisation rates as well as spa- 
tial deviations in the IGM temperature due to the inho- 
mogeneity of the cosmic web supp ort this view at least on 
scales of few co - moving Mpc (see iMaselli fc Ferraral 120051 : 
lFurlanettdl2009l : iMeiksml 12009] and references therein). On 
larger scales (~ 100ft -1 Mpc co-moving) the source spa- 
tial distribution and their spectral variabilit y could be a n 
additional cause of variations in the UVB fZuol Il992bl lal: 
iMeiksin and Whitdl2003t iBolton fc Vie1l201lft . 

The fluctuations induced by radiative transfer effects 
could also be efficiently recorded in the ionisation state of the 
metals, because their rich electronic structure and atomic 
spectrum are more sensitive, compared to H and He, to 
the radiation field fluctuations l|Ohll2002l : IFurlanetto fc Loebl 
120031 : lFurlanettdl2009T ) . 

Numerical schemes which solve the cosmological radia- 
tive transfer equation by applying differe nt approximations 
are n ow quite mature and well tested (|lliev et al. 2006a. 
2009 and references therein for an overview of the avail- 
able codes) and are able to simulate complex scenarios in- 
volving 



Ciardi et 



large cosmolo gical boxes and number of source s (e.g . 
it al.ll2003alh~];lniev et all2006d;lTr"ac and Cen Il2007i: 



McQuinn etai1l200d : lBaek et al.N20ld : ICiardi et all 12012ft 



Typically, these codes are restricted to the hydrogen chem- 
istry, with only a few of them including a self-consistent 
treatment of the helium component, which is particularly 
relevant for a correct deter mination of the gas temperature 
(see e.g. ICiardi et alll2012ft . None of them though includes 
the treatment of metal species. 

In the interstellar medium (ISM) community, on the 

other hand, severa l photo-ionisati on codes, as Cloudy 

jFerland et alJll99Sft. MAPPINGSI II (| Allen et al.lll99Sft and 
MOCASSIM (|Ercolano et al.l 12003ft are able to simulate the 
complex physics of galactic HII regions largely polluted by 
heavy elements. 

The aim of the present work is to des cribe a novel exten- 
sion of the radiativ e trans fer code CRASH (| Ciardi et al 1 l200ll : 
IMaselli et aT]|2003l . lioogj ). which has been integrated with 
Cloudy, allowing the prediction of metal ionisation states 
self-consistently with the radiation field as calculated by the 
radiative transfer through the IGM density field. 

The paper is structured as follows. In Section 2 and 3 
we briefly introduce the radiative transfer code CRASH and 
the photo-ionisation code Cloudy. In Section 4 we illustrate 
the details of the integration of the two codes into a self- 
consistent pipeline. The tests of the new CRASH variant called 
CRASH3 are reported in Section 5. Section 6 summarizes the 
conclusions. 



2 CRASH 

CRASH is a 3D radiative transfer (RT) code designed to fol- 
low the propagation of hydrogen ionising photons (i.e. with 
energy E 13.6 eV) through a gas composed by H and He. 
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The code adopts a combination of ray tracing and Monte 
Carlo (MC) sampling techniques to propagate photon pack- 
ets through an arbitrary gas distribution mapped on a carte- 
sian grid, and to follow in each grid cell the evolution of gas 
ionisation and temperature. This treatment guarantees a re- 
liable description of such evolution in a large variety of con- 
figurations, as shown by the Cosmological Rad iative Trans- 
fer Comparison Project tests (|lliev et al.ll2006ah and its var- 
ious applications to th e stud y of the H and He reionisation 
jCiardi et alj l2003al lbl. l2012h . the imprints of the fluctuat - 
ing background on the Lya f orest |Maselli fc Ferrarall2005h . 
the quasar proximity effects l|Maselli et al.1 12007, 2009]) and 
the impact of reionisation on the visibility o f Lya emitters 
|Daval et al.ll201ll ; [jeeson-Daniel et al.ll2012T l. 

The MC algorithm adopted allows to easily add new 
physical processes. In its first version (jCiardi et all l200ll ) 
the code describes H photo-ionisation due to point sources, 
and includes the effect of re-emission following gas recom- 
bination. The subsequent versions brought about signifi- 
cant improvements. First the physics of He and the thermal 
evolution of the gas have been introduced, together with 
the treatment of an ionising b ackground field (jMaselli et al.l 
20031: iMaselrife Ferrarall2005l ). In the latest release (CRASH2; 
Masem et aT1l2009l . hereafter MCK09) we introduced multi- 
frequency photon packets obtaining significant improve- 
ments in terms of accuracy of the ionisation and temper- 
ature profiles, as well as computational speed. Hereafter the 
code name CRASH will refer to the version CRASH2. 

In parallel with the reference code, variants and 
extensions have been developed such as: MCLya 
l|Verhamme et all [2006), which has adapted the refer- 
ence algorithm to treat the resonant propagation of Lya 
photons; CRASHa IjPierloni et al.l 120091. which follows the 
self-consistent propagatio n of both Lya phot ons and ionis- 
ing continuum radiation; IPartl et al. I l|201ll ) have instead 
developed an MPI parallel implementation of the base 
CRASH algorithm. 

The new release described in this paper (hereafter 
CRASH3) extends the standard CRASH photo-ionisation algo- 
rithm to the treatment of the most cosmologically relevant 
metals in atomic form: C, O and Si. The current numeri- 
cal scheme is based on a new pipeline which combines the 
excellent capabilities of CRASH in tracing the radiation field 
with the sophi sticated features of the photo-ionisation soft- 
ware Cloudy l|Ferland et al.l ll998). The inclusion of a con- 
siderably large set of data imposed by the numerous metal 
ionisation states has required a substantial source code re- 
engineering that introduces a new memory management 
and a re-modelling of the photon packet-to-cell interaction. 
CRASH3 is consequently more modular, optimized and easily 
integrable with other codes. In the following, we briefly re- 
view the basic ingredients of the CRASH algorithm that are 
required to understand the implementation of the new ver- 
sion, and we refer the reader to the original papers for more 
specific details. 

A CRASH simulation is defined by assigning the initial 
conditions (ICs) on a regular three-dimensional cartesian 
grid of a given dimension cells) and a physical box linear 
size of Lb, specifying: 

• the number density of H (tih) and He (riHe), the 
gas temperature (T) and the ionisation fractions (xhii = 



riHIl/riH, XHell = WHell/MHe and a?HeIII = TlHelll/ftHe) at the 

initial time to; 

• the number of ionising point sources (N 3 ), their posi- 
tion in cartesian coordinates, luminosity (L s in ergs -1 ) and 
spectral energy distribution (SED, S s in erg s _1 Hz -1 ) as- 
signed as an array whose element provide the relative inten- 
sity of the radiation in the correspondent frequency bin (see 
MCK09 for more details); 

• the simulation duration tf and a given set of intermedi- 
ate times tj £ {to, . . . ,tf} to store the values of the relevant 
physical quantities; 

• the intensity and spectral energy distribution of a back- 
ground radiation, if present. 

The simulation run consists in emitting photon packets from 
the ionising sources and following their propagation through 
the domain. Each photon packet keeps traveling and deposit- 
ing ionising photons in the crossed cells, as far as its content 
in photons is completely extinguished or it escapes from the 
simulated box (although periodic boundary conditions in the 
packets propagation are possible). 

At each cell crossing, CRASH simulates the radiation-to- 
gas interaction by evaluating the absorption probability for 
a single photon packet as: 



P(t) = 1 -e 



(1) 



where r is the total gas optical depth of the cell given by 
the sum of the contributions from the different species, i.e. 
t — thi + thoI + THcii- The number of photons absorbed in 
the cell can then be estimated as: 

iV 7 (1 - e~ T ) , (2) 

where A r 7 indicates the photon content of a packet entering 
the cell. 

The number of the deposited photons for each spectral 
frequency is then used to compute the contribution of photo- 
ionisation and photo-heating to the evolution of xmi, iHeii, 
SHein and T. The s et of equations solv ed in the code is 
described in detail in lMaselli et al.l (|2003l ) and in MCK09. 

For the sake of the following discussion, we remind here 
that the evolution of the thermal state of the gas in each 
cell is regulated by the formula: 



dT 
~dt 



3fcsp 



k B T^+U(T, Xi)-A{T,Xi 



(3) 



where p is the number of free particles per unit volume, H 
and A are respectively the heating and the cooling func- 
tions and the subscript i refers to all the ion species com- 
posing the gas. Ub is the Boltzmann constant. The heating 
function is determined by computing the photo-heating re- 
sulting from the photon-to-gas interaction discussed above, 
while A is calculated by adding up the contribution of vari- 
ous radiative processes: collisional ionisation and excitation, 
recombinations, Bremsstrahlung and Compton cooling. Dif- 
ferently from photo-heating, these processes are treated as 
continuous process es, described by their respective rates (see 
iMaselli et al1l2003l for details). 



3 CLOUDY 



Cloudy (|Ferland et al.|[l998l ) is a code designed to simulate 
the physics of the photo-ionised regions produced by a wide 
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class of sources ranging from the high temperature blue stars 
to the strong X-ray emitting Active Galactic Nuclei. The 
main goal of Cloudy is the prediction of the physical state of 
photo-ionised clouds including all the observably accessible 
spectral lines. The latest stable release of Cloudy (at the 
time of writing v. 10. cQ) simulates a gas which includes all 
the heavy elements of the typical solar composition and the 
contribution of dust grains and molecules present in the ISM. 

In this Section we will focus on the description of the 
Cloudy features that have been primarily used to implement 
CRASH3. The reader interested in the details of the code im- 
plementation or in reviewing the many physic al processes in- 
cluded can find mor e appr opriate references in lFerland et al.l 
|l998l ) and lFerlandl (|2003l ). 

Unlike CRASH, Cloudy is a ID code assuming as pre- 
ferred geometrical configuration a symmetrical gas distri- 
bution around a single emitting source, with photons prop- 
agating along the radial direction. Cloudy can also simu- 
late the diffuse continuum re-emitted by recombining gas 
as nearly isotropic component under the assumption that 
the diffuse field contribution is generally small and can be 
treated by lower order approximations. Additional isotropic 
background fields can also be handled, as long as their shape 
and intensity are specified by the user. Some popular back- 
ground models (like the Haardt and Madau cosmic UV spec- 
trum described in iMadau fc Haardtll2009r) are already dis- 
tributed with the code. The contribution of the Cosmic Mi- 
crowave Background (CMB) radiation can also be accounted 
for because it is an important source of Compton cooling for 
low density gas configurations typical of the IGM. 

The micro-physics implemented in Cloudy is very ac- 
curate: it include s all the metals present in t he typical so- 
lar composition (|Grevesse and Sauval I ll99St) described as 
multi-level systems and treated self-consistently with the 
ions of the lightest 30 elements. Photo-ionisation from va- 
lence, inner shells and many excited states, as well as col- 
lisional ionisation by both thermal and supra-thermal elec- 
trons and charge transfer, are included as ionisation mecha- 
nisms. The gas recombination physics is simulated including 
the charge exchange, radiative recombination, and dielec- 
tronic recombination processes. Cloudy simulates all these 
processes adopting an approximation method for the radia- 
tion field evalu ation known as escape probabilities method 
(|Hubenvl200ll ), instead of evaluating the full radiative trans- 
fer as done by CRASH. This choice implies the loss of many 
details pertaining the line properties description, e.g. line 
pumping by the incident continuum, photon destruction by 
collisional deactivation and line overlap. In the standard re- 
lease of CRASH the details of the lines are not accounted for 
and the above limitations are thus negligible. 

Once the characteristics of the source and the species 
involved in the calculation are set up, Cloudy estimates 
the radial distribution of the ionisation and temperature 
fields by simultaneously solvi ng the equations of ionisa- 
tion and thermal equilibriu m l|Dopita fc Sutherland! [20031 ; 
lOsterbrock fc F crland 2005|). A static solution describing the 
physical state of the gas is then found by dividing the do- 
main in smaller regions and iterating until convergence is 
reached. The usual assumption of such calculations is that 
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atomic processes occur on timescales much shorter than the 
temporal scales of the macro-physical system. Cloudy does 
not 'a priori' assume that the gas is in equilibrium and 
the solution provided is a general non-equilibrium ionisation 
configuration for a static medium that does not retain any 
information of the temporal evolution of the system towards 
the converging state. A large set of information about the 
relevance of the physical processes that establish the final 
convergence, and the details of the line emission processes 
are also provided with a great level of detail. This micro- 
physic treatment can not be directly handled in CRASH with 
the same level of accuracy. 

Cloudy describes the thermal equilibrium of the photo- 
ionised gas providing the local thermal balance obtained in 
each simulated sub-region. In the absence of non-thermal 
electrons produced by high-energy photons, this thermody- 
namic equilibrium is generally specified by the electron tem- 
perature T e , because the electron velocity distribution of the 
gas is predominantly Maxwellian. Despite CRASH does not 
distinguish between the contribution of the different species 
to the gas temperature as calculated in eq. [31 in practice 
T ~ T e to a good approximation, and thus it is consistent 
with the output from Cloudy. It should be clarified though 
that Cloudy provides a gradient of temperatures within the 
simulated region. Because the CRASH resolution is a single 
cell, whenever we use a temperature provided by Cloudy 
this should be intended as the volume average over a cell. In 
the following we will always refer to the temperature as T. 

Implementing CRASH3 we have minimised the differences 
between the two codes, e.g. by disabling in Cloudy all the 
physical processes non explicitly treated in CRASH such as 
molecule and dust grain physics, charge transfer effects and 
radiation pressure, among others. In addition, while CRASH 
simulates the propagation of hydrogen ionising photons, 
Cloudy requires that any spectral information is provided 
in the energy range 13.6136 • 10" 8 eV< E < 100 MeV. 
For this reason, the spectrum used as input for Cloudy is 
the one provided by CRASH in the frequency range 13.6 eV^ 
E ^ Emax (see following Section), while it is set to zero for 
13.6136 • 10" 8 eV< E <13.6 eV and E max < E <100 MeV, 
where E max is the energy corresponding to the higher fre- 
quency sampled in the CRASH simulation (see MCK09 for 
more details) . The deeply different set-up of the simulations 
in CRASH and Cloudy, both in geometrical configuration and 
in the gas micro-physics, has made the building of the com- 
mon pipeline a non trivial task. This will be detailed in the 
following Sections. 



4 THE INCLUSION OF METALS IN CRASH 

The extension of the CRASH algorithm with the full micro- 
physics of metals is hardly viable because of the extreme 
complexity of the metal electronic structure which would 
increase the computatio nal costs of a RT simulation t o unac- 
cepted) ly high levels fsee lDopita fc Suth erland 2003; Drainc 
1201 ll and references therein for a modern introduction to the 
physics of the metals in ionised regions). 

For this reason, we have adopted a hybrid approach in 
which CRASH performs the RT only through H and He, while 
the metal ionisation states are implemented self-consistently, 
but they are calculated with Cloudy. The two algorithms 
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interact within a single pipeline called CRASH3, which is 
sketched in Figure[T] It should be noted that, if the radiative 
transfer is performed in a cosmological context, the pipeline 
applies to each single redshift, z. In this case, in addition 
to the ICs of the RT simulations, other physical quantities 
might depend on z, e.g. the cooling off the CMB radiation. 

More specifically, CRASH3 recognises the metal enriched 
sub-domain by applying a masking technique to it and by 
labelling the enriched cells, which are indicated here with 
the subscript m. After masking the m-cells containing met- 
als, it derives the spectral energy distribution and luminos- 
ity of the ionising radiation present in each of the m-cells 
(S,L) m , which are then used to query a dynamic database 
(DB) of configurations pre-computed with Cloudy to iden- 
tify the corresponding ionisation states of the metals and 
the temperature of the gas. This procedure is repeated for 
a set of times tk £ {to, ...,tf} and the physical state the 
enriched cells is then evaluated as a sequence of kf equi- 
librium configurations at tk- To ensure that the ionisation 
level of the metals is consistent with the energy configura- 
tion, the pipeline checks that the ionisation fractions of H 
and He evaluated by Cloudy and CRASH are in agreement. 
An example of this internal convergence test can be found 
in Section (2H7TJ 

It should be noted that this approach neglects the con- 
tribution of metals to the optical depth, as well as to the 
diffuse radiation emitted by recombining gas. This approxi- 
mation is justified as long as the metal abundances are very 
small compared to those of H and He. In fact, when the 
number density of heavy elements nz ~ 10 _3 nH, the total 
photo-ionisation cross section starts to be dominated by the 
metal component fsee lDraindl201l1 . Chapter thirteen for a 
complete discussion). Because the metallicity observed in 
the IGM is typically below this value (see the Introduction), 
the above assumption is justified in the cases of our interest. 

As detailed in Section 3, Cloudy can deal with either 
a single source configuration or a background radiation. If 
only a single source is present in the computational domain 
of CRASH or when each cell is illuminated from a preferential 
direction, then the Cloudy point source configuration should 
be used to estimate the ionisation and temperature state 
of the cell. When instead a cell is illuminated more or less 
uniformly from all directions, then the background radiation 
set-up should be adopted. 

In the current implementation of CRASH3 we have in- 
cluded the species C, O and Si, which are the most abun- 
dant metals observed in the IGM (see the Introduction). 
However, the inclusion of other species, which may be rel- 
evant for a more accurate estimate of the gas cooling func- 
tion, is a straightforward extension of the present scheme 
(see Sec. 5.1.5). 

Despite the conceptual simplicity of the approach de- 
scribed above, the coupling of CRASH and Cloudy in a sin- 
gle numerical scheme presents several technical challenges. 
In the following we will describe in more detail the strat- 
egy adopted for such integration and some aspects of the 
pipeline. 



4.1 Step 1: initial conditions 

The starting point of the pipeline is the set-up of the CRASH3 
ICs, which are the same of CRASH as specified in Section 2, 
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Figure 1. CRASH3 simulation pipeline. The quantities S, L, Xi 
and T are the SED, luminosity, ionisation fractions (per species 
i) and gas temperature. The subscripts c and m refer to all the 
cells in the computational domain and the metal enriched cells 
only, while the superscripts c and C3 refer respectively to the 
quantities provided by the radiative transfer simulation in Step 
2 and by the database (DB), i.e. as a product of the full CRASH3 
pipeline. See text for more details. 



with the addition of the spatial distribution and abundance 
of all the metal species present in the computational domain. 
These can be artificially created by hand or can be obtained 
as a result of e.g. hydrodynamic simulations that include 
physical prescriptions for metal production and spreading. 

A preliminary analysis of the spatial distribution of met- 
als allows to identify the m-cells that need to track the ra- 
diation field (Step 2) and then query the pre-computed DB 
(Step 3). To identify these cells a boolean mask is built iso- 
lating the enriched portion of the simulation volume. The 
building of the mask can be performed before the begin- 
ning of the simulation and passed as additional IC or can 
be created in memory during the simulation initialization. 
If the mask contains m-true values, a shadow map of m cell 
spectra, S£, and luminosities, is allocated to store the 
shapes and luminosities of the incoming packets: each time 
a packet enters a cell, the mask is used to check whether the 
cell is an enriched one and consequently and L„ should 
be updated (see Step 2). The masking has been designed to 
optimise the performances of the code, but in principle the 
SED and luminosity can be calculated in each cell of the 
computational domain, if required. 



4.2 Step 2: radiative transfer simulation 

The next step (Step 2 in Fig. [1]) consists in performing a RT 
simulation which, in addition to the evaluation of (xi,T)^ in 
all the cells of the domain, tracks also the SED and luminos- 
ity of the ionising radiation in each of the m metal enriched 
cells (5*, L)m- All the above physical quantities are stored at 
times tk, as already mentioned above. The values of (S, L)^ 
at tk are determined adding up the contribution of all the 
incoming photon packets up to that time. In practice, the 
code tracks all the multi-frequency packets entering the cell 
and calculates as the sum of the luminosities of all the 
packets, while Sm as the average SED in each frequency bin. 
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4.3 Step 3: determination of the physical state of 
the metal enriched cells 

Finally, a search engin^fl looks for the Cloudy precom- 
puted configuration that best matches the values of H, He, 
metal number densities and (xhii, ZHeii, :EHeiii, S, L)^(tfc). 
The best fitting configuration then provides the ionisation 
fractions of the metal ions and the temperatures T^ 3 of the 
gas in the metal enriched sub-domain (Step 3 in Fig. [TJ. 

If the matching criteria are not satisfied (see below and 
Appendix B for more details) , the database is extended with 
additional on-the-fly Cloudy runs. It is important to point 
out that Step 3 does not severely affect the basic algorithm 
performances because it is confined only to the m enriched 
cells and a large number of Cloudy calculations are pre- 
computed and stored in the DB. The interested reader can 
find more details about the lookup strategy in Appendix B. 

A correct computation of the temperature is not trivial, 
because also in the absence of metals the temperatures pre- 
dicted by CRASH and Cloudy are not in perfect agreement in 
the vicinity of the point sources (MCK09). More generally, 
it has been shown that different approaches to the radia- 
tive transfer do not always predict consistent temperatures 
in such regions (|lliev et al. I l2006bh . This means that every 
time a discrepancy between and T^ 3 occurs it is impor- 
tant to understand if this is due to the presence of metals or 
just to the differences in the two codes. For these reasons we 
allow some flexibility in the matching criteria for the tem- 
perature and we define a temperature deviation AT^ as: 

AT' = T c — T cs (4) 

where i = met refers to the deviation calculated for a gas 
contaminated by metals, while i = pris refers to a pristine 
gas. The difference AT m = AT™ et - AT£" S is > by design 
and it is due only to the effect of metals. In the enriched cells 
in which AT m is greater than some threshold value for the 
maximum tolerated deviation, is replaced by T„ 3 . We 
apply this correction (accounting for the metals feedback on 
the temperatures) when the AT m is higher than ~10% T^ 3 
without metals. However the value adopted for the threshold 
depends on the physical problem being study and for this 
reason the criteria for the threshold can be defined at the 
beginning of each simulation. Note that the temperature 
correction has some weak feedback on the physical state of 
the gas also via its recombination coefficients, especially for 
the helium component. 

In realistic applications, the temperature provided by 
the hydrodynamics (Thydro) could be significantly higher 
than the temperature established by the RT. This gener- 
ally happens in gas shocked regions where the metal ions 
are determined mainly by collisional ionisation rather than 
photo-ionisation. Collisional ionisation is correctly handled 
by the pipeline at Step 2 for H and He, because CRASH se- 
lects T = max(T c , Thydro) during its simulation initialisa- 
tion process, maintaining the right temperature. Note that 
these regions are known as ICs, and a masking technique can 



2 We have implemented a Standard Query Language database 
using one of the codes freely available to the scientific community. 
More details can be found in Appendix B 



be used to isolate them from the standard pipeline work- 
flow. The ionisation status of the metal component eventu- 
ally present in these cells can be determined from Cloudy 
pre-computed tables, and by interpolating the gas number 
density, metallicity, temperature, as well as the ionisation 
fractions of hy drogen and helium. For mor e details on this 
technique see lOppenheimer fc Pavel (2006) and references 
therein. 



5 TESTS 

In this Section we present three tests designed to establish 
the reliability of the new code. The first test (Section I5.1[> 
concentrates on the standard Stromgren sphere case, albeit 
of a metal enriched gas. The second test (Section l5.2|) in- 
vestigates the sensitivity of CRASH3 to fluctuations of the 
radiation field induced by different source types and tracked 
by the many metal ionisation states. Finally, the third test 
(Section I5.3[) describes a more realistic physical configura- 
tion by using as ICs those from a snapshot of a hydrody- 
namic simulation. 

Hereafter the gas metallicity (or equivalently the metal 
mass fraction) is defined as Zg — Mz/M g , where Mz is the 
total mass of the elements with atomic number higher than 2 
and M g is the total mass of the gas; the metal mass fraction 
in the Sun is set to Zq w 0.0126, according to the metal 
abundances relative to hydrogen as reported in the Cloudy 
Hazy guide I and references therein, and taking into account 
the 10 most abundant elements: H, He, C, N, O, Ne, Si, Mg, 
S and Fe. 

Unless otherwise stated, the metal feedback on the cal- 
culation of the gas temperature is neglected and we use 
k f = 100. 



5.1 Test 1: Stromgren sphere with metals 

We consider a configuration similar to the one in Test 2 
of the Cosmologic al Radiative Transfer Comparison Project 
(llliev et al. 2006a), but for the presence of metals. 

The simulation box has a co-moving side length of 
6.6 kpc and it is mapped onto a regular grid of iV 3 = 128 3 
cells. The gas is assumed to be uniform and neutral at 
the initial temperature T = 100 K, with a number den- 
sity rigas = 0.1 cm -3 and a hydrogen (helium) number 
fraction of 0.9 (0.1). Only one point source is considered 
with coordinates (1,1,1), a black-body spectrum at temper- 
ature T = 10 5 K and ionisation rate of N = 10 51 phot g- 1 
(i.e. a luminosity L ~ 5 • 10 40 ergs -1 ). To ensure a good 
convergence of the MC scheme, the source radiation field 
has been sampled by 2 • 10 s photon packets. The redshift 
of the simulation is set at z = and the simulation dura- 
tion at tf = 5 • 10 8 yrs, i.e. several times the recombination 
time characteristic of the simulated gas configuration. This 
choice assures that convergence is reached at the end of the 
simulation. We uniformly contaminate the gas with carbon 
(n c ~ 2.2 ■ 10" 7 cm" 3 ), oxygen (no ^ 4.41 ■ 10~ 7 cm" 3 ) 
and silicon (nsj — 3.12 ■ 10~ 8 cm -3 ), corresponding to 
Z g = 6.45 • W~ 3 Z e . These numbers are obtained maintain- 
ing the metal abundances relative to hydrogen mentioned 
above. Finally, the source spectrum is sampled by 91 fre- 
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Figure 2. Convergence between Step 2 (dashed lines and vari- 
ables with superscript c ) and Step 3 (solid lines and variables 
with superscript < - 73 ) for Test 1 in the absence of metals (see text 
for details). The lines refer to the simulation time tf = 5 ■ 10 s yrs. 
At distances larger than 70 cells the gas is neutral and there- 
fore it is not shown in the plots. Top panel: profile of xhi 
(green lines) and rrjjn (blue). Middle panel: same as above for 
^Hel (red lines), xjjeii (gray) and XHelll (black). Bottom panel: 
Axi = — x? 3 for j=Hll (blue line), Hell (gray) and He III 
(black). 

quencies and it extends to an energy E max — 0.2 keV, in 
order to include the photo-ionisation edge of the ion O vi. 

Although we contaminate the entire box with metals, it 
is possible to significantly reduce the number of DB queries 
taking advantage of the spherical geometry of the problem 
and assuming that each radial direction is equivalent. This 
is justified as long as the number of photon packets used 
is large enough that the fluctuations induced by the Monte 
Carlo sampling are negligible. We then apply a spherical av- 
erage of all the relevant physical quantities at each distance 
d (expressed in cell unit) from the ionising source, and we 
use these averaged values to search for the better-matching 
configuration in the DB. All the quantities discussed and dis- 
played in this Test should be intended as explained above, 
i.e. as spherically averaged values. 

5.1.1 Pipeline convergence 

Before proceeding further with the analysis of the results, 
we discuss the internal convergence of the pipeline. More 
specifically, we need to assure that the ionisation fractions 
a^Hii, ZHeii and iHein calculated in Step 2 and Step 3 of the 
pipeline are in agreement with some user-defined threshold. 
The same needs to be true for the temperature T if the 
effect of metals is negligible. This would ensure that the 
physical configuration and energetic balance evaluated in 
Steps 2 and 3 are fully consistent. To this aim we run Test 1 
in the absence of heavy elements. 

In Figure [2] we show the profile of :ehi, £hii (top panel), 
KHelj KHeii and CEHelll (middle panel) as evaluated by Step 
2 (dashed lines and variables with superscript c ) and Step 
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Figure 3. Top panel: Temperature convergence between Step 

2 (T c , dashed line) and Step 3 (T C3 , solid line) for Test 1 in 
the absence of metals (see text for details). The lines refer to 
the simulation time tf = 5 ■ 10 s yrs. Bottom panel: relative 
difference AT /T = (T C3 - T C )/T C3 . 

3 (solid lines and variables with superscript ° 3 ) at the time 
tf = 5 • 10 8 yrs. The values of xm and xmi result in agree- 
ment within 10 -4 for d < 40 cells (~ 2 kpc) and within 10~ 2 
up to the ionisation front (I-front), identified as the location 
where the ionised fraction drops below ~ 0.8. The agreement 
degrades to ~ 18 % (with respect to the maximum absolute 
difference between ionisation fractions) in the two cells in 
which the curves of xm and xmi cross, and then it restores 
up to 10~ 3 . Both codes predict the crossing in the same cell. 
A similar behaviour is shown in the middle panel for he- 
lium. A discrepancy of ~ 7 % in the XHell curves is seen at 
a distance d ~ 22 cells, when the abundance of Hen starts 
to increase. He in results in an even better agreement, up to 
10 -4 . Reasonable accuracy (less than 20 %) is also reached 
in the fronts of Hen and Hei where both algorithms pre- 
dict a similar shape. Because at the end of the HII region 
the intensity of the radiation has dropped substantially (see 
following Section), sometimes Cloudy does not find a conver- 
gent solution in few iterations, and this results in a poorer 
agreement between Step 2 and 3. 

The above differences are more evident in the bottom 
panel of the Figure, where we show the absolute difference 
Axi = xf — xf' 3 for i=Hn, Hen and He in. 

Despite the satisfying agreement in the CRASH3 imple- 
mentation, some small discrepancies remain due to the dif- 
ferences between the CRASH and Cloudy geometries and the 
numerical implementation of the ionisation and energy equa- 
tions (see Sec. 3). We found that a critical ingredient to reach 
an acceptable convergence is to sample the source spectrum 
with a large number of frequency bins (see details in the Test 
Section). This is necessary because the helium component is 
very sensitive to this sampling, in particular in the vicinity 
of the ionisation potential of He n. 

The temperature radial profiles estimated in the two 
Steps are shown in Figure [3] (top panel). The curves present 
large discrepancies in the cells near the source, with a rel- 
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Figure 4. Normalized logarithm of the spherically averaged spec- 
tra. Lines refer to the spectra at a distance d (from the top to the 
bottom) of 20, 27, 38, 45, 55, 69 and 70 cells. The spectra corre- 
spond to a time tt = 5 • 10 8 yrs. As a reference, some ionisation 
potentials are also indicated as vertical lines. 



ative variation AT/T = (T cz - T C )/T C3 as high as 60 
percent. The difference drops below 30 percent for d > 30 
cells. This is not reflected in the H and He profiles shown 
in Figure [2] because of the weak temperature dependence 
in the gas recombination coefficients. Such a difference has 
been already noticed and discussed in the CRASH vs Cloudy 
comparison test in MCK09 and can be ascribed just to the 
different implementation of the temperature estimate in the 
two codes. 

Because CRASH updates the temperature (compared to 
its initial value) only in those cells reached by ionising pho- 
tons, outside the HII region T drops to the initial value of 
100 K. On the other hand, the temperature calculated by 
Step 3 is provided by Cloudy and, as already mentioned 
above, all the regions where the intensity of the radiation is 
too faint do not provide a convergent solution. In the few 
cells in which Cloudy cannot provide a reliable estimate be- 
cause of the too weak radiation field, the CRASH temperature 
T c is retained as better estimate of the physical tempera- 
ture at the front. 



These convergence tests have been repeated using dif- 
ferent ICs for the gas number density and the source lu- 
minosity. More specifically, we have run simulations on a 
grid of cases with values n gas = 1, 0.1, 0.01 cm -3 and 
N = 5 ■ 10 50 , 5 ■ 10 51 phot s -1 . It is found that, as the gas 
density decreases, the agreement improves for the H species, 
while small discrepancies still remaining in the He species. 
Such discrepancies are, on the other hand, always below 20 
% and remain limited to the small number of cells encom- 
passing the ionisation fronts. A similar improvement is ob- 
served also for the gas temperature convergence. 

Hereafter all the variables will refer to values calculated 
at Step 3 and the superscript Ci will be omitted to simplify 
the notation. 



E io „ [eV] 
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Si 


E x i 


13.598 


24.587 


11.260 


13.618 


8.152 


Exii 




54.400 


24.383 


35.118 


16.346 


Exin 






47.888 


54.936 


33.493 


E x iv 






64.494 


77.414 


45.142 


E x v 






392.090 


113.900 


166.770 


E x vi 






489.997 


138.121 


205.060 



Table 1. Ionisation potentials for H, He, C, O and Si up to the 
ionisation level VI as used in Cloudy. 
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Figure 5. Logarithm of the photo-ionising luminosity, Li on , as 
function of d at a time tf = 5 ■ 10 s yrs. Different curves are calcu- 
lated integrating all the photons above the ionisation threshold of 
(from the top to the bottom curve): Hi (black line), Hel (blue), 
Sim (pink), Siiv (gray), Cm (brown), Hen (violet), Civ (red), 
OlV (green), Ov (orange), Ovi (yellow). 

5.1.2 Ionisation field 

In Figure [4] we show how the spectral shape of the ionising 
radiation field described in terms of its spherical average 
(obtained as described at the beginning of Sec. 5) changes 
with the distance d as a result of geometrical dilution and 
filtering. Each line refers to the simulation time tf = 5 ■ 
10 8 yrs. 

The spectra shown are truncated at E — 120 eV to 
have a better visualisation of the most relevant line poten- 
tials. The upper curve corresponds to a distance of d — 20 
cells, while the lower to d — 70 cells; at larger distances the 
intensity of the radiation becomes too faint to solve the ion- 
isation equilibrium in every configuration at Step 3 of the 
CRASH3 pipeline. Because the medium is fully transparent in 
H and He up to a distance d = 20 cells, we have normalized 
the spectra in the Figure to the maximum value of the spec- 
trum at d = 20 cells. This choice allows a better visualisation 
of the shaping effects by the partially ionised regions. The 
ionisation potentials of the metals enriching the box (see Ta- 
ble [TJ are also shown as reference even if, by design, metals 
do not contribute to the filtering of the ionising radiation. 
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Figure 6. Fractions of the various components as a function of distance d from the source in Test 1. The values are taken at the 
simulation time tj = 5- 10 8 yrs. From the top panel to the bottom the species are: H (green lines) and He (blue); C (red); O (orange); and 
Si (black). In each panel the same ionisation states are represented by the same line-styles: solid lines refer to the neutral components 
(e.g. Ol), long dashed to the first ionisation state (e.g. Oil), short dashed to the second (e.g. Oin), dotted to the third (e.g. Oiv), long 
dashed-dotted to the fourth (e.g. Ov), short dashed-dotted to Ovi, and dashed-spaced to Ovu. 



The absorption due to Hi, Hei and Hen is clearly visible 
in correspondence of the respective ionisation potential, i.e. 
13.6 eV, 24.6 eV and 54.4 eV. 

Figure[5]shows the total photo- ionising luminosity, Lion, 
available for ionisation of some reference species as a func- 
tion of d at a time tf = 5 • 10 8 yrs. Li OTl is defined as: 

-Bmax 

Lion (d) = h^ 1 J S{E,d)dE, (5) 

where Ei on is the ionisation potential of the species consid- 
ered (see Table [T] for reference) and h p is the Planck con- 
stant. 

Because the spectrum includes only photons with E > 
13.6 eV, Lion is an underestimate for those elements with 
an ionisation potential below 13.6 eV, i.e. carbon and sili- 
con. Note that Loi, Lcn and Loin overlap respectively with 
Lm (black line), Lnei (blue) and LhcII (violet) because of 
the very similar ionisation potential. If we increased the fre- 
quency resolution of the spectrum the curves would show 
some small difference. This would be at the expense of the 
computational time without significant advantages in the 
accuracy of the metal ionisation state. For this reason we do 
not further increase the frequency resolution. Notice that all 
the luminosities decrease smoothly with d; this is consistent 
with the behaviour expected in an HII region, as already re- 



ported and extensively commented in iMaselli et al.l (|2003l ) 
and MCK09 for the hydrogen and helium components. 



5.1.3 Metal ionisation states 

We now analyse the behaviour of the metal ionisation states, 
plotted in Figure [6] Before discussing the details of the Fig- 
ure, it is necessary to point out that the balance among the 
different ions is primarily established by the relative values 
of their ionis ation potentia ls and of their recombination co- 
efficients (see lFerlandll2003T ) . It also depends on the spectral 
distribution of the radiation field and its variations with the 
distance from the source, as induced by the radiative trans- 
fer effects. The complex interplay between these numerous 
processes makes the interpretation of the results non trivial; 
despite this, some trends have a straightforward explana- 
tion. 

Because of the small amount of metals included in this 
test, we expect their impact on the evolution of H and He 
to be negligible. This is confirmed by a comparison of the 
curves in the top panel of the Figure to the correspond- 
ing curves in Figure [2] which are obtained without metals. 
The maximum difference is ~ 7 % across the I-front of Hn. 
Additional effects on the ionisation fractions induced by an 
increase in the gas metallicity will be investigated in Section 

E33] 
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Figure 7. Fractions of the various components as a function of the distance d from the source in Test 1. From the top panel to the 
bottom the species are: H and He; C; O; and Si. In each panel colors refer to different simulation time: t = 10 6 (blue), 10 7 (red) and 
5 ■ 10 8 yrs (black). 



We now turn to analyse the behaviour of carbon (sec- 
ond panel from the top). For d < 30 cells, C is in the form of 
Cm, Civ and Cv, with a predominance of the latter close 
to the source. This high ionisation level is obtained from 
a combination of collisional ionisation and photo-ionisation. 
The evolution of C ill is very similar to that of He n because 
of the similar ionisation potentials (see Fig. 3}. The abun- 
dance of Civ is dictated by the evolution of Cm and Cv, 
and their relative recombination coefficients. C iv is present 
throughout the entire HII region, with xciv being always 
below 30 percent. For d > 30 cells a;cv goes to zero because 
only few C iv ionising photons are available (see Fig. [5] as 
a reference). The ionisation potential of Cv is outside our 
frequency range (see Table [TJ and thus higher ionisation 
states are not present. Because of the paucity of photons 
with E > -Ecm, at d > 30 cells only Cm is present in 
large quantities with xcm ~ 70%. At d ~ 60 cells, similarly 
to what happens to Hn and Hen, also Cm declines and 
C ii dominates. Finally, outside the HII region, only C i is 
present. 

In the third panel from the top the ions of the oxygen 
are shown. The ionisation potential for Ovi is the highest 
photo-ionising energy available in the adopted spectrum. A 
very small fraction of O vn is in fact present in few cells 
around the source. It should be noticed that collisional ioni- 
sation contributes to this fraction for ~ 10 % at T ^ 7T0 4 K, 
which is present when d < 3 cells. The presence of Ovi, is 
more evident but it is consistently limited to the inner re- 



gion of the ionised sphere and decreases rapidly with the 
distance from the source in favour of lower ionisation levels. 
For d > 30 cells Oiii dominates the ionisation balance un- 
til it drops in favour of O n, roughly at the end of the HII 
region. As for the other species, outside the HII region only 
Oi is present. 

In the bottom panel the behaviour of silicon is reported. 
Siv completely dominates the inner region of the Stromgren 
sphere, with a long tail extending to d ~ 55 cells, where it 
is in equilibrium with lower ionisation states. In the central 
region many ions are in equilibrium with a low ionisation 
fraction. Sim reaches a maximum fraction of ~ 0.4 at the 
center of the HII region. Sin dominates at d > 52 cells. The 
abundance of Sii in the outskirts of the HII region is much 
lower than the one of e.g. Ci because of the lower number 
density of Si and the much larger cross section of Sii at the 
resonance (see Fig. 13.2 of lDraindlamh . 

In general, it is possible to say that in the vicinity of the 
source the most abundant species are those with the higher 
ionisation state compatible with the maximum potential in 
the spectrum, i.e. H n, He in, C v, O v and Siv. Despite Eoy, 
Eoyi and i?siv being covered by the spectrum, the abun- 
dance of photons at these energies is so low that xovi, zovn 
and xsivi are negligible. As the distance increases, the lumi- 
nosity available for ionisation decreases, in particular for ions 
with high ionisation potentials (see Fig. This is reflected 
by the decrease of the abundance of these highly ionisation 
states and the predominance of lower ionisation states (e.g. 
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Figure 8. Fractions of C (upper panel), O (middle) and Si (lower) ions as a function of the distance d from the source in Test 1. 
In each panel, lines with crosses refer to Test 1 with the source spectrum extending to a lower ionising energy of Esn, rather than Em. 
Note that the neutral fractions in the former case are equal to zero in the range of cells shown here. 



Hen, Cm, Om and Sim). Species like Siiv and Civ are al- 
ways present, although they are not dominant, because the 
spectrum of the ionising radiation maintains energies higher 
than -Esim and Ecm throughout the HII region (see Fig.JS]). 
At even larger distances (d > 60 cells) the dominant species 
are typically singly ionised metals and the neutral compo- 
nents, due to the absence of ionizing radiation. 

While the discussion above refers to the final gas con- 
figuration, in Figure [7] we show some reference species at 
different simulation times, i.e. t = 10 6 , 10 7 and 5- 10 8 yrs. It 
should be kept in mind that in this case no metal feedback is 
included in the calculation (see Sec. 15.1.5]) and thus the pro- 
files of the metal species at each time are independent from 
each other. While the shape of each species changes sub- 
stantially between 10 6 and 10 7 yrs, for t £ 10 7 yrs they have 
almost converged. This is especially true for the metal ioni- 
sation potentials larger than the third and for He in. Overall, 
the comments made for Figure [H] at t/ = 5 ■ 10 8 yrs apply 
also at earlier times. 



5.1-4 Effect of photons with E < 13.6 eV 

In this Section we discuss the impact of photons with E < 
13.6 eV in the evaluation of metal ionisation states. Because 
the ionisation potential of Ci and Sii is below Em, i.e. in a 
range which is not covered by CRASH and where the spectrum 
is set by default to zero, this means that photo-ionisation 
by photons with E < 13.6 eV is neglected, resulting in a 



systematic underestimate of xcu and xsm as visible in the 
outer region of the Stromgren sphere. 

To investigate this limitation of our pipeline, we have 
extended the spectrum to a lower energy of Esn = 8.152 eV. 
It should be kept in mind, though, that in the range Esu-Em 
the luminosity of the spectrum is overestimated because the 
absorption by metals (the only species contributing to the 
optical depth in this energy range) is not accounted for. This 
means that the change in the spectrum at these frequencies 
is due only to geometrical dilution. The approximation is 
more severe beyond the Hen front (see Figure [5}, where 
the only surviving radiation has energy below Em and the 
optical depth of the medium is dominated by the metal^j. 

In Figure [5] we summarise the results of this run by 
showing the evolution of different C, O and Si ions. Lines 
with crosses refer to the case in which the source spectrum 
extends below 13.6 eV. In the inner part of the HII region the 
agreement between the curves with and without the inclu- 
sion of the low energy tail is excellent for the neutral frac- 
tions and the first ionisation state, while it gets worse for 
higher ionisation states. The reason is that, when the spec- 
trum is extended below 13.6 eV additional physical processes 
become relevant, in particular line emission from collisional 



3 In realistic configurations other absorbers like dust and 
molecules could play a dominant role in establishing the opti- 
cal depth of the medium at these frequencies, but these are not 
accounted for in CRASH3. 



© 2010 RAS, MNRAS 000, ??42T1 



12 L. Graziani et al. 



excitation. The most prominent case is the collisional exci- 
tation of the O iv line. While the behaviour of iov is the 
same in both cases (and for this reason it is not reported 
here), in the presence of non Hi- ionising radiation Oiv re- 
combines more quickly to O ill, inducing the discrepancies 
observed between ~ 20-35 cells. Something similar happens 
to C v and C III, while C iv is not affected because it is not 
a dominant species at these distances, as shown in Figure [6] 
For the Si the differences are smaller. The picture in the out- 
skirts of the Hll region instead changes completely: not only 
do xcn, xou and aisiii extend outside the Hll region, while 
the corresponding neutral components disappear, but also 
the higher ionisation states are affected. This results from a 
combination of the photo-ionisation due to the photons with 
energies below 13.6 eV and the additional physical processes 
mentioned above. At these distances the former effect has a 
larger impact than in the vicinity of the source, where ionisa- 
tion is dominated by photons with energies above 13.6 eV. 
It should be noted that outside the HII regions, or, more 
correctly, when the absorption of photons is not dominated 
any more by H and He, our assumption of neglecting the 
metal contribution to the optical depth fails and the metal 
ionisation states are not calculated correctly any morfl 



5.1.5 Feedback by metals 

In this Section we investigate the effect of metals on the 
evaluation of the gas temperature. In the Appendix A the 
convergence of this feedback with respect to the choice of 
kf is investigated. This is done by calculating T at the time 
tf = 5 • 10 s yrs in the standard configuration of Test 1 and 
then changing the metallicity of the gas. 

Test 1 has been repeated changing the gas metallicity 
Z g , while maintaining the relative abundance of C, O and 
Si. In the upper panel of Figure [9] we show the results in 
terms of: 



AT/T 



T(0)-T(Z g 
T(0) 



(6) 



where T (0) is the value of the temperature relative to a con- 
figuration with Z g — 0. The reference case (Z g = 0.0064^0, 
dashed line) does not show any significant metal cooling, 
with the exception of the region near to the source (d < 5 
cells) where recombination and re-emission of high ionisation 
states of C, O and Si is more significant, inducing an average 
AT/T ~ 10%. Temperature deviations at the Hen I-front 
(d > 65 cells), where the ionising radiation is very faint, are 
also present, with AT/T < 10%. Increasing the metallicity 
to one percent solar (solid line) does not change the results. 
Only at Z g — 0.064Z© (dashed-dotted line) some cooling is 
visible at each distance. In few cells near the source AT/T is 
as high as 40%, while it remains below 10 percent at d > 10 
cells. The effect of metal cooling starts to be very significant 
for Zg ^ 0.64Z© (dashed-spaced line), with a AT/T ~ 60% 
in the inner part of the HII region and ~ 15 — 20% also in 
its outer part. 



4 Other important processes playing a relevant role in the out- 
skirts of the HII region, as e.g. + H < — > O + , have been 
deactivated in Cloudy by disabling the charge t ransfe r. More ref- 
erences on the subject can be found in iDrainei J201lh . 
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Figure 9. AT/T = (T(0) - T(Z 9 ))/T(0) as function of the 
distance d from the source in Test 1 at the simulation time 
tf = 5 ■ 10 s yrs (see text for more details). Upper panel: the 
curves refer to a gas enriched by C, O and Si, and a metallicity 
Zg of: 0.0064Zq (dashed line, reference value), O.OIZq (solid), 
O.O64Z (dashed-dotted), O.lZg (dotted) and O.64Z (dashed- 
spaced). Lower panel: the curves refer to a gas enriched only 
with C, O and Si (black lines) or with the ten most abundant 
elements in the solar composition (red). For both cases the gas 
metallicity Z g is O.OIZq (solid lines) and O.IZq (dotted). 




30 40 
d [cells] 

Figure 10. Fraction of Civ as function of the distance d from 
the source in Test 1 at the simulation time tj = 5 • 10 s yrs. The 
curves refer to a gas enriched only with C, O and Si (black lines) 
or with the ten most abundant elements in the solar composition 
(red). For both cases the gas metallicity Z g is O.OIZq (solid lines) 
and 0.1Z Q (dotted). 



Even if it is not shown in Figure [9] we have computed 
the solar and super-solar cases (Z g ^ IZq), finding a devi- 
ation of AT/T > 50% in the inner region, and values ex- 
ceeding 40 percent up to d ~ 45 cells. These results confirm 
the dominant role played by the metal cooling in this metal- 
licity range. These cases though should be considered only 
as indicative because in this metallicity regime the metal 
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Figure 11. Sketch of the geometrical set-up used for Test 2. 



contribution to the absorption cannot be neglected and the 
assumptions of our method fail. 

As a second step we have investigated the dependence of 
our results on the gas composition; this is shown in the bot- 
tom panel of Figure^] More specifically, while the black lines 
are the equivalent of the ones in the upper panel, the red 
lines have been obtained including the 10 most abundant el- 
ements in the solar composition. Note that the red and black 
lines are not distinguishable in the reference metallicity case. 
In both cases the number densities of metals are such that 
their abundance relative to H is the same as those in the 
solar composition. It is clear that, at a fixed metallicity, 
the contribution of C, O and Si to the cooling is dominant 
compared to other elements, such as N, Ne, Mg, S and Fe. 
We have thus neglected such elements in all further tests. 
It should be noted that just adding the contribution of N, 
Ne, Mg, S and Fe to the gas metallicity (without keeping 
it constant) would bring Z g from e.g. 0.0064^0 to O.OIZ©, 
and O.O64Z to 0.lZ e . 

In Figure [10] we finally report the effect of a varying 
metallicity and chemical composition on the ionisation frac- 
tion of Civ. It is clear that it is safe to neglect metals 
other than C, O and Si as their impact on the global ionisa- 
tion equilibrium is limited at a maximum of few percent for 
Z g ^ O.IZq, while it is close to zero for lower metallicity. A 
similar effect is observed for all the other relevant ions. On 
the other hand, the dependence on the gas metallicity Z g is 
higher and it changes for the various species. In some cases, 
like the one shown in the Figure (as well as e.g. O in and 
Siiv), the ionisation fraction of the species increases with 
increasing metallicity, while in others (as e.g. C in, O v and 
Siiv) the trend is reversed. It is beyond the scope of this 
paper though to discuss this issue in more details. 



5.2 Test 2: metal fluctuations in the overlap of 
two HII regions 

Here we study the behaviour of metals in the HII region over- 
lap produced by two point sources located in cells (20, 64, 64) 
and (40,64,64). This geometrical set-up is sketched in Fig- 
ure [TT] and it is designed to investigate the sensitivity of 
our method to small changes in the source characteristics. 
This is done by varying the source ionising rates, N\ and 
N2, and their associated black-body spectral temperatures, 
Ti and T2. All the other numbers are the same as those in 
Test 1, with the exception of the gas number density which 
is n gas — 1 cm -3 , to obtain a sharper overlap profile. Be- 
cause in this test we are interested in studying the behaviour 
of metals only in the overlap region, we concentrate on the 
plane corresponding to x — 30 (blue circle in Fig. Ill[) . 




h [cells] 



Figure 12. Ionisation fractions as function of h, evaluated in 
the plane equidistant from the two sources (x = 30) of Test 2 in 
the reference case. From top to bottom the panels refer to: iehii 
(green dashed) and £Hell (blue dashed); xrjn ( re d dashed), iciii 
(red solid), 23111 (black dashed), £sim (black solid), xrjn (orange 
dashed) and xrjin (orange solid); xqiv (red dotted), ajgnv (black 
dotted) and rrg;v (black dashed-dotted). 



5.2.1 Reference case 

We first describe the set-up used as reference case, with two 
identical sources with N\ — N2 — 9 ■ 10 51 phot s _1 and Tj = 
T2 = 10 5 K. Similarly to Test 1, here we show the results 
by averaging the input physical quantities on all the cells of 
the plane at the same distance h from cell (30,64,64) (see 
Fig. HH> . The results refer to the final time tf = 5 ■ 10 8 yrs. 
Figure[l2]shows the profiles of the ionization fraction of some 
selected ions. 

The profile of the H and He species is very similar to 
that of a single Stromgren sphere (see Test 1), with the 
exception of ^Heiii, which is always confined in regions very 
close to the sources and thus is not shown. 

In the middle panel Cn, Cm, Sin, Sim, On and Oin 
are shown together because they trace the external regions 
of the two overlapping Stromgren spheres (see Fig. [6]). Cm, 
Oin and Sim are present for h < 12 cells with different 
ionisation fractions (zoiii ~ 1, icm ~ 0.5 and xsuu ~ 0.3), 
indicating that these ions have a different sensitivity to the 
ionising field; this is also in qualitative agreement with the 
relative trends noticed in Test 1 (see Fig. [5]). 

In the lower panel a; civ, ^Siiv and xsiv are shown. 
Similarly to the species analysed in the previous panel, 
they present an almost constant value throughout the fully 
ionised region, with the exception of xsiv, which starts de- 
clining earlier, in favour of lower ionisation states. The ab- 
sence of O iv, O v and C v (which reach an ionisation fraction 
of only a few percent) is consistent with the absence of He ill. 
Compared to the corresponding species in Test 1, here the 
ionisation fractions are higher due to the larger ionisation 
rate. 
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Figure 13. Ionisation fractions as function of h, evaluated in the 
plane equidistant from the two sources (x = 30) of Test 2 when 
the source ionisation rate is changed compared to the reference 
case (see text for details). From the top to the bottom the panels 
refer to mmi, ^ciVi a^silll an d :ESiiv- In all panels solid lines refer 
to the reference case with N\ = N 2 = 9 ■ 10 51 phot s — 1 ; dashed to 
TVx = 9 ■ 10 51 phot s" 1 and N 2 = 7 ■ 10 51 phot s -1 ; short-dashed 
to Ni = 9 ■ 10 51 phot s" 1 and N 2 = 3 ■ 10 51 phot s _1 ; dashed- 
dotted to N± = N 2 = 7 ■ 10 51 phot s -1 ; dotted to JVi = N 2 = 
3 ■ 10 51 phot s" 1 . 



5.2.2 Variations in the source ionisation rates 

In this Section we discuss the variations in the results in- 
duced by changes in the ionisation rates. These are shown 
in Figure [13] for some reference species. 

First we have decreased the ionisation rates of the two 
sources simultaneously, maintaining the symmetry of the 
problem, i.e. using Ni = N 2 — 7 ■ 10 51 phot s _1 and 
3- 10 51 phot s _1 , the latter being the minimum ionising rate 
allowing for an overlap of the two HII regions. We have then 
decreased the ionisation rate of only one source, while main- 
taining iVi = 9T0 51 phot s _1 . As expected, the extent of the 
fully ionised region (upper panel) decreases with decreasing 
total ionisation rate. A similar trend is observed also in the 
profile of the other species, which changes smoothly with de- 
creasing ionisation rate. Compared to the behaviour of xhii 
though, these differ in that also the amount of ionisation de- 
creases. Remarkably, a;siiv seems to be insensitive to small 
variations of the ionisation rates, at least in the inner parts 
of the HII region. This suggests that other species would be 
more suitable to give indications on the characteristics of 
the sources. 

By comparing qualitatively the results with those in 
Figure [6] we can conclude that CRASH3 produces the pre- 
dictable behaviour for the ions also when more than one 
source is present and that it is sensitive to small changes 
in the source ionisation rates. This is very important for 
realistic applications of the code. 



Figure 14. Ionisation fractions as function of h, evaluated in the 
plane equidistant from the two sources (x = 30) of Test 2 when 
the source spectrum is changed compared to the reference case 
(see text for details) . From the top to the bottom the panels refer 
to xhii, ^CIVj ^Silll and I'snv- In all panels solid lines refer to 
the reference case with Ti = T 2 = 10 5 K; the dashed to Ti = 
T 2 = 5 • 10 4 K and the short-dashed to Ti = T 2 = 5 ■ 10 5 K. 



5.2.3 Variations in the source spectra 

In this Section we study the variations induced by changes in 
the temperatures Ti and T 2 of the black-body spectra. The 
source ionisation rates are set to the reference value Ai = 
A2 = 9- 10 51 phot s _1 . The results in terms of distribution of 
some reference ionisation fractions are shown in Figure [141 

In the top panel the H 11 I-front is identical in the cases 
Ti = T 2 = 5 ■ 10 4 K and Ti = T 2 = 10 5 K, while it recedes 
significantly for T\ = T 2 — 5 ■ 10 K. This is because by 
increasing the temperature from 5 • 10 4 K to 10 s K, most 
of the photons still have E < i?Heii and get preferentially 
absorbed by hydrogen because of its higher abundance. For 
these two cases iHeii (which is not plotted here) has a profile 
very similar to that of xuu, with the only difference that it 
slightly increases with increasing temperature because more 
photons have E > -EhcI- On the other hand, when the tem- 
perature of the black-body spectra is increased to 5 • 10 5 K, 
most of the photons have an energy in the vicinity of -BhcII, 
causing a drop in zhii and iHeii and the appearance of XHein 
(which is not plotted here), which was missing in the other 
cases, being relegated only in the immediate surroundings 
of the sources. 

While for Ti = T 2 = 5 ■ 10 4 K there is hardly any 
C iv available because of the lack of enough photons with 
E > -Ecin, as the spectrum temperature increases so does 
xciv- Increasing the black-body temperature further, and 
thus shifting the peak of the black-body spectrum to higher 
frequencies, brings a drop in xciv, which is balanced by an 
equal increment of lev (which is not shown in this plot). 

Unlike in the previous tests, in this case both Si in and 
Si iv show remarkable changes by varying the spectrum tem- 
perature. While zsiiii keeps decreasing as the temperature 
increases and higher ionisation states are favoured, xsav in- 
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Figure 15. Slice cut through the simulation box showing the 
initial distribution of log(nci) [cm -3 ] in Test 3. 

creases from Ti = T 2 = 5 • 10 4 K to T\ = T 2 = 10 5 K, while 
it decreases when T\ — T2 — 5 ■ 10° K, similarly to xciv- 

Also in this case we can conclude that CRASH3 is very 
sensitive to changes in the spectrum of the ionising sources. 

We have verified that a similar conclusion applies when 
only one of the two source temperature is changed or when 
we have adopted a power-law rather than a black-body spec- 
trum. 



5.3 Test 3: radiative transfer on a cosmological 
density field enriched by metals 

In this Section we present an extension of the Test 4 pro- 
posed in the Cosmologica l Radiative Transfer Comparison 
Project (|lliev et al . 2006a). The original test set-up has been 
adopted, but we have extended the ICs to include He, C, O 
and Si. For reference, the box size is L b = 0.5/i _1 Mpc co- 
moving, h — 0.72, N c = 128, the simulation duration is set 
to tf = 4 • 10 s yrs, starting at redshift z = 9. The cosmo- 
logical evolution of the box is disabled and the simulation 
starts with a neutral gas at initial temperature T — 100 K. 

The H and He fractions are the same as in Test 1. 16 
point sources are present in the box with a black-body spec- 
trum at T = 10 5 K. The sampling of the spectrum is done 
with 91 frequencies, to reach the accuracy established in our 
Test 1. Finally, we have used 10 s packets to ensure a good 
convergence as done in MCK09. It is important to note that 
the aim of this test is to show how the CRASH3 pipeline is 
applied to a realistic density configuration and the rich set 
of information provided by it, and it is not meant to quan- 
titatively reproduce observational results. 

Differently from the previous tests in which we popu- 
lated the full computational volume with metals, here we 
assign a metallicity of Z g — A • 0.0064Zq to those cells with 
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Figure 16. Slice cut through the simulation box showing the dis- 
tribution of the gas temperature T in the same slice of Figure 1151 
at the time t = 5 ■ 10 4 yrs in Test 3. 

an overdensity A = p/{p) > 10. This results in about 5 per- 
cent of metal enriched cells. The number density of C, O and 
Si relative to the hydrogen are the same ones introduced at 
the beginning of the Test Section. 

In Figure [15] the distribution of the C i number den- 
sity is shown in a slice of the simulation box containing the 
brightest source, located within the densest region in the 
computational volume. The distribution of Si i and O i is 
the same, with nsa = 0.142 ■ nci and noi = 2 • nrji. By con- 
struction, most of the metals are concentrated in the vicinity 
of the sources and in general in the highest density regions. 

As shown in Figure [9] of Test 1, gas with metallicity 
Z g > O.IZq starts to cool with different efficiency depend- 
ing on the distance from the illuminating source (see com- 
ments in the Section), while for Z g > 0.64Zq the cooling 
is relevant at every distance from the source. In the box of 
this Test, cells with A > 100, are only 0.012 percent of the 
total number, while cells with A > 50 are about 4 percent. 
We then expect the effect of metal cooling to be negligible 
on the global ionisation status of the medium and we ignore 
it in discussing the results. In realistic 3D simulations the 
metal cooling should be carefully accounted for in all the 
cells with Z g > 0.1Zo, because it acts differently depending 
on the gas metallicity and the RT effects shaping the field. 
The statistical effects of the metal cooling will be analysed 
with care in future applications in the context of physically 
motivated enrichment patterns. 

To better understand the distribution of the metal ion- 
isation species discussed in the following, in Figure [16] we 
show the distribution of gas temperature in the same slice 
as above at a time t = 5 • 10 4 yrs. Black regions (T ~ 100 K) 
correspond to areas not reached by the radiation field, the 
green color clearly traces the He 11 ionisation fronts at a typ- 
ical temperature T ~ 2 • 10 4 K, while the red/orange is as- 
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Figure 17. Slice cut through the simulation box showing the distribution of (from left to right and from top to bottom): rrcVi ^CIV; 
^Cllli 2;oV; 2;oiVi ^Ollli ^SiVi ^SilV an d ^Silll- The slice shown is the same of the previous Figures at the time t = 5 ■ 10 4 yrs in Test 3. 



sociated to regions at T ~ 4 ■ 10 4 K, dominated by Hn and 
Hen. Finally, warmer, yellow regions correspond to areas 
where there is also substantial Hem. 

The distribution of some reference metal species is 
shown in Figure [T7] across the same slice. While in the vicin- 
ity of the brightest source C v is present in large quantities, 
with lev as high as 1, only small traces are visible further 
away. On the other hand, these regions are typically rich 
in C in, showing a qualitative behaviour similar to the one 
discussed in Test 1, where Cv and Cm are complementary 
species. Consistently with the same Test 1, Civ is present 
everywhere along the filaments with relatively low abun- 
dances. Regions with fractions close to zero are those which 
have not been reached by ionising radiation (see Fig. I16[) . 
A similar behaviour is observed for O, with Ov present in 
very small traces only in the vicinity of the brightest source, 
O iv extending a bit further away along the filaments and 
O in being present in substantial amounts elsewhere. Finally, 
while there are only traces of Siiv, a;siv is as high as 1 in the 
vicinity of the source and Si in is present in modest quanti- 
ties everywhere. 



While it is not possible to make a quantitative compari- 
son with Test 1, the qualitative results in terms of ionisation 
fractions at a given distance from the brightest source are in 
good agreement, with the highest ionisation levels present 
only in the close proximity of the source and the lower ion- 
isation levels present everywhere along the filaments hit by 
the ionising radiation. 

Similarly to what done in Section r5.1.4l we have investi- 
gated the impact on the evaluation of metal ionisation states 
of photons with E < 13.6 eV by extending the spectrum to 
to energies lower that 13.6 eV down to -Esu- In Figure fTSl the 
differences in ionisation fractions (Aa;; = Xi t E sn — Xi,s H i) 
computed in the two set-ups are reported for the same ions 
of Figure [171 As for Test 1, we find that where xmi ~ 1 
the agreement between the case with and without the in- 
clusion of the low energy tail is excellent in most of the 
metal enriched domain, xorv and zoiii present an average 
discrepancy of about 13%, with a maximum of ~ 40% con- 
fined in few cells where the ionisation fractions are relatively 
low (i.e. below ~ 30%). For higher x the agreement is excel- 
lent, reflecting the results in Section [57L4] where the largest 
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Figure 18. Slice cut through the simulation box showing the distribution of the differences in ionisation fractions obtained in Test 3 
with the reference source spectrum and with a spectrum extending to a lower ionising energy of i?sii (Aii = Xi,.E g .j — £i,E HI ; see text 
for more details). From left to right and from top to bottom the panels refer to the species: Cv, Civ, Cm, Ov, OlV, Oiii, Siv, SilV 
and Si III. The slice shown is the same of the previous Figures at the time t = 5 ■ 10 4 yrs. 



discrepancies were observed in correspondence of those cells 
in which the abundance of the different species was rais- 
ing or declining. On the other hand iov always shows a 
good agreement between the two cases. The Si and C ioni- 
sation fractions present a maximum discrepancy below 10% 
for their higher (i.e. iv and v) ionisation species, while it 
raises to ~40% for zciii and xsan- Also for these species 
only few cells have such large discrepancies. 

In general we can conclude that, similarly to what seen 
in Test 1, the extension of the spectrum to energies lower 
than 13.6 eV results in some discrepancies in the evalua- 
tion of the metal ionisation fractions. Such discrepancies are 
limited to some of the metal ions (especially the lowest ion- 
isation states) and are largely due to a number of physical 
processes which are neglected (see footnote of Test 1 and 
the discussion in Section [5.1.4[) . In the selected slice, they 
are also confined to a small number of cells typically located 
far from the ionising sources. As a reference, the percentage 
of metal enriched cells which present a discrepancy larger 



than 5%, 10%, 25% and 50% in the visualised species is 
19.3%, 9.6%, 4.2% and 0.6%, respectively. 



6 CONCLUSIONS 

In this paper we presented CRASH3, the latest release of the 
3D radiative transfer code CRASH. In its current implemen- 
tation CRASH3 integrates into the reference algorithm the 
code Cloudy to evaluate the ionisation states of metals, self- 
consistently with the radiative transfer through the most 
abundant species, i.e. H and He. The feedback of the heavy 
elements on the calculation of the gas temperature is also 
taken into account, making of CRASH3 the first 3D code for 
cosmological applications which treats self-consistently the 
radiative transfer through an inhomogeneous distribution of 
metal enriched gas with an arbitrary number of point sources 
and/or a background radiation. It should be noted that the 
algorithm is based on the assumption that metals do not 
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contribute to the gas optical depth, which is valid in typical 
configurations of cosmological interest, as e.g. in the IGM. 
On the other hand, the code is not suitable for applications 
in which metallicities £ Z@ are involved. 

The pipeline has been tested in idealized configurations, 
such as the classic Stromgren sphere, albeit enriched with 
heavy elements, as well as in a more realistic case of multiple 
sources embedded in a polluted cosmic web, which shows 
the rich set of information provided by the metal ionisation 
states. Through these validation tests the new method has 
been proved to be numerically stable and convergent with 
respect to a number of parameters. 

The dependence of the results on e.g. the source char- 
acteristics (spectral range and shape, intensity), the metal 
composition, the gas number density and metallicity has 
been investigated. Despite the difficulty in testing quanti- 
tatively the results, we find that the qualitative behaviour 
is consistent with our expectations. 

In conclusion, CRASH3 is an excellent code to simulate 
the evolution of various species in a low metallicity gas il- 
luminated by an ionising radiation, which can be modelled 
both as a background radiation and point sources. CRASH3 
will then be an invaluable tool to investigate the interpre- 
tation of e.g. metal absorption lines in quasars' spectra or 
fluctuations in the UVB, in more detail, with a better mod- 
elling of the relevant physical processes and with a higher 
accuracy than done before. These applications are topics for 
further investigations. 
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APPENDIX A: CONVERGENCE TESTS 

When the temperature feedback is introduced in the CRASH3 
pipeline (see description of Step 3) kf times, the RT process 
is altered. While the gas recombination depends weakly on 
the temperature and we can expect negligible effects, the nu- 
merical convergence of the metal ionisation fractions needs 
be verified for different values of kf. Here we show the results 
of some tests run to check the convergence of our results with 
respect to kf when the metal feedback is enabled. 

In Figure IA1I we show the profile of the variation of 
a reference ionisation fraction (xcu) and temperature for 
Test 1 run with the metal feedback switched on. The lines 
refer to the final time tf = 5 ■ 10 8 yrs and to values of 
kf — 10, 50, 100 and 500. The convergence of the results 



© 2010 RAS, MNRAS 000, ??-f2T1 



20 L. Graziani et al. 



within the ionised sphere (d ^ 40 cells) is excellent: even for 
small values of kf the relative discrepancies remain within 
few percent both in the ionisation fraction (top panel) and 
in the temperature (bottom panel). Numerical oscillations of 
a few percent (or smaller) in Ax/x are present in a handful 
of isolated cells, but they disappear for kf ^ 100. 

Closer to the I-front (d ^ 50), where :rcn is very small 
(i.e. sC 10 -4 ), the convergence is worse, in particular in the 
2-3 cells encompassing the I-front. This behaviour was ex- 
pected because, as already discussed in Test 1, the accu- 
racy of our scheme is diminished in the outer regions of the 
ionised bubble, where the Cloudy computation is not always 
convergent. It should be noticed, however, that xcn presents 
one of the worse behaviours and higher ionisation states are 
more numerically stable. We have checked that the conver- 
gence is similar also at earlier times, when the ionised sphere 
has not reached an equilibrium configuration. Other species 
show a similar dependence on kf. For example, xciv presents 
an excellent convergence (of the order of a percent) within 
the ionised region already with kf = 10, while in the prox- 
imity of the I-front, where xciv ^ 10 -6 , it degrades to tens 
of percent. 

The temperature shows an even better behaviour, with 
a relative convergence of a few percent even in the vicinity 
of the I-front. Also at earlier times the convergence is never 
worse than 10%. The oscillations in the curves have the same 
origin of those in the upper panel. 

We should notice that the value of kf needed to reach 
a good convergence is likely dependent on the problem at 
hand. In general, though, we can conclude that we expect 
a value of kf ~ 500 to be sufficient for most applications 
which require the feedback of metals on the determination 
of the gas temperature. 



APPENDIX B: DATABASE IMPLEMENTATION 
AND LOOKUP STRATEGY 

Here we provide more details on the implementation of the 
relational database introduced in Section 4 and the lookup 
strategy we use to query its data. 

The introduction of a relational database storing pre- 
computed data is motivated by the huge number of Cloudy 
computations (N comp ) needed by our pipeline in specific 
metal configurations. Consider a general reionisation sim- 
ulation involving N z density snapshots in which the metals 
are computed tf times on a sub-domain of N m (z) cells; the 
total number of computations (N cornp ) is: 

N comp = N z x Nm(z) xt f . (Bl) 

Note that N comp depends indirectly on the adopted grid 
resolution N c and the metal filling factor, therefore the size 
of the metal map N m can rapidly reach a high value. For 
instance, a simulation with N c = 128 3 , N m ~ 10 5 (i.e. cell 
filling factor of 5%) and tf — 10 requires N cornp ~ 10 6 com- 
putations in a single snapshot. As the direct evaluation of 
N cornp Cloudy runs during the RT simulation is feasible just 
in few cases, we use a DB of precomputed configurations to 
alleviate the total computational cost. 

The DB stores a large number of Cloudy runs 
abstracted as object-data associating initial conditions 



(RUNICs) with the resulting ionisation fraction and tem- 
perature (RUNRESULTS); this association is uniquely 
identified by an indexing key (RUN_ID). 

The set of RUN_ICs is composed by the number densi- 
ties of all the species n m polluting the metal contaminated 
sub-domain, the luminosity L m and the spectral shape S m , 
specified on a fixed number of frequency bins (see Section 
4). As detailed in the Test Section, 91 bins are generally re- 
quired to obtain the accuracy of the tests in matching the 
HII fronts. 

The RUNRESULTS set is composed by the ionisation 
fractions of all the species x m :i and the temperature T CA . 

As prescribed by the relation database theory, all the 
variables are organised in columns of the database tables 
and related by the RUN ID. The database solution is prac- 
tically implemented by combining two relational databases: 
a reduced "in- memory DB" (called production DB) used 
during the RT (see Step 3 of the Pipeline in Section 4.3), 
and an "archiving DB" (a standard TCP/IP accessible DB) 
collecting the results of many runs and increasing in time. 
Before running a simulation the content of the smaller DB 
is extracted from the archiving DB and customised on the 
problem at hand as explained below. 

The content of the production DB can be tuned on 
the specific problem under investigation because the CRASH3 
pipeline is applied in post-processing and part of the 
RUN_ICs can be partially pre-constrained. First, the n m 
values are exactly known from the hydro-simulation at fixed 
redshift, while the values of L m and S m can be only pre- 
dicted because they are affected by the metal feedback. Sec- 
ond, in the regime of applicability of our pipeline we tested 
that the metal cooling feedback is weak on the RT field and 
by running a RT simulation just through H and He (the only 
absorbers, see Section 4) we can have a reasonable estimate 
of the statistics of the produced spectral shapes. Note also 
that in absence of feedback from metals L m and S m are pro- 
vided exactly from Step 2 at every time but N cornp can be 
huge and still require some pre-computed configurations to 
successfully finish the run. 

The DB is then populated by spanning a subset of com- 
patible ICs providing a non-redundant set of final ionisation 
fractions and temperatures; this is realised by testing the 
sensitivity of the photo-ionisation software to small varia- 
tions (Sn, SL, SS) of their values. An example of these tests 
is reported in Section 5.1 when we studied the effects of the 
metal feedback on the resulting temperature and ionisation 
fractions (see Figures 9 and 10). 

During a CRASH3 simulation, the pipeline is respon- 
sible for lookup the DB at Step 3 as explained in Section 
4.3. To recognise a compatible pre-computed configuration 
in each m-cell of the metal polluted sub-domain, Step 3 per- 
forms a cascade of SQL-SELECT queries by comparing the 
RUN ICs of the interested cell with the RUN ICs records 
stored in the DB tables, as explained below. 

The first query (Q_l) selects a ResultSet (R_l) by 
querying for the compatible number densities: n m ± 5n m . 
The number of records (i.e. table rows) in R_l is defined as 
Ri. This is realised by a series of SQL-BETWEEN conditions 
on n m , concatenated by .AND. clauses; in SQL meta-code: 

SELECT FROM < TABLE NAME > 
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WHERE ( 

(nHI BETWEEN nHI + - dnHI) 
AND 

(nHel BETWEEN nHel + - dnHel) 
AND . . . 

) 

A second query Q_2 is performed on R_l applying 
the same logic to the cell luminosity L m and providing a 
ResultSet R_2 (R 2 < Ri) in which L m G L m ± 5L m . 

Finally, a Q_3 query finds compatible spectral shapes 
S m in the ResulSet R_2 by simultaneously comparing the 
shape in every spectral bin within the tuned oscillation 8S. 
Q_3 provides a ResultSet R_3 with R3 <C 7?2- This last 
query involves all the spectral bins and drastically reduces 
the number R3 of final candidate configurations. On the 
other hand it is necessary to strictly compare the spectral 
shapes in every bin to guarantee that the spectral fluctu- 
ations due to RT effects are correctly accounted for and 
the full compatibility between results coming from the DB 
queries and the ionisation fractions produced by the com- 
plementary on-the-fiy runs. 

The set R_3 of candidate pre-computations is further 
reduced to R3 = 2 final candidates by minimising the dif- 
ference of their H and He ionisation fractions with the val- 
ues found at Step 2 (See Section 4.3). We generally allow a 
maximum difference below 20 % for H and 10 % for both He 
fractions, but these values can be adapted to the problem at 
hand, increasing the precision of the matching criterion. The 
final result is then found as mid-point of the linear interpo- 
lation of their H,He and metal fractions if they are either 
sides of the Step 2 results, otherwise just the configuration 
with the closer value is considered, first checking for arHeii 
and XHein, and finally for xhii- As an example of the ac- 
curacy obtained in the resulting configuration consider the 
convergence test in Section 5.1.1 for a HII region. 

Also notice that in all the cells with xmi = ^Heii = 
£Hein = 1 (e.g. the first 8 cells of Figure 2), the configu- 
rations in the set R 3 cannot be disentangled. In this case 
either the metal ionisation fractions differ within a threshold 
value (normally 10 -3 ) and one of them is accepted, or the 
entire R_3 set is rejected and a on-the-fly run is performed. 

As final comment we point out that the DB solution has 
the only purpose of reducing the computational complexity 
of the problem by reusing previous computations and cannot 
provide any general solution of the RT through metals, but 
must be targeted on the specific problem at hand. Also note 
that for a maximum accuracy in the metal ion predictions, 
and when the value of N comp implies a reasonable computa- 
tional time, the DB matching thresholds can be reduced or 
the DB lookup can be excluded from the pipeline allowing 
just on-the-fiy computations. This is the case of most sin- 
gle snapshot computations with medium resolution and low 
metal filling factor. 
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